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ABSTRACT 

Using an Eulerian perturbative calculation, we show that the distribution 
of relative pairwise velocities which arises from gravitational instability of 
Gaussian density fluctuations has asymmetric (skewed) exponential tails. The 
negative skewness is induced by the negative mean streaming velocity of pairs 
(the infall prevails over expansion), while the exponential tails arise because 
the relative pairwise velocity is a number, not volume weighted statistic. The 
derived probability distribution is compared with N-body simulations and shown 
to provide a reasonable fit. 

Subject headings: large-scale structure of universe — galaxies: interactions 
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1. Introduction 

Redshift surveys present a distorted picture of the world because peculiar motions 
displace galaxies from their true spatial positions. This phenomenon, which would make 
redshift surveys useless for intergalactic spaceship navigators, is extremely useful for 
cosmologists. It can serve as a probe of the dynamics of gravitational clustering and 
the cosmological mass density parameter, Q ( [Sargent and Turner 19771 ; Peebles 1980| , 



hereafter LSS; [Kaiser 1987| ; [Hamilton 1992| ; [Peebles 1993| , hereafter PPC; [Regos and Szalay 



1995). A convenient statistical measure of the distortion effect is the galaxy two-point 
correlation function in redshift space. Under certain assumptions it can be expressed as 
a convolution of the true spatial correlation function, £(r), with the distribution of the 
relative line-of-sight velocities of pairs of galaxies, p(w\r,9) . Here r and w are respectively, 
the spatial separation and relative radial velocity of a pair of galaxies, while 9 is the angle 



between the separation vector r and observer's line of sight (cf. LSS; |Fisher 1995| , hereafter 



F95). The purpose of this Letter is to derive p(w\r, 9), using weakly nonlinear gravitational 
instability theory. This distribution was measured from N-body simulations and estimated 
indirectly from redshift surveys. At r ;$ 1 /i _1 Mpc, whereQ the galaxies are strongly 
clustered (£ ^ 20), the observations are consistent with an exponential distribution ( [Peebles 



1976] |Davis and Peebles 198$ |Fisher at al. 1994| , hereafter F94; [Marzke et al. 1995| ; |Landy. 



Szalay, fc Broadhurst 1997] ) . The fact that p(w) at small separations differs strongly from 
its initial, Gaussian character, is not surprising: after all, the small-scale velocity field has 
been 'processed' by strongly-nonlinear dynamics in clusters, and exponential distributions 



were recently derived from the Prcss-Schcchtcr (1974)| theory ( Sheth 1996 . Diafcrio and 



Gcllcr 1996| ). On larger scales, where the fluctuations have small amplitudes, one naively 



expects to see the 'unprocessed' initial conditions. However, N-body experiments suggest 



3 We use the standard parametrisation for the Hubble constant, H = 100 h km s Mpc . 
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that p(w\r, 6) retains its exponential character even at separations r £ 10 /i _1 Mpc, where 
£ £0.1, despite the fact that the initial density and velocity fields in those experiments were 
drawn from a Gaussian distribution (Efstathiou et al. 1988| , hereafter EFWD; [Zurek et~aT 



1994| , hereafter ZQSW; F94). At similar separations, an exponential p(w\r, 6) has also been 



inferred from observations (F94; |Loveday et al. 1996|) . The simulations also show that the 
radial component of the distribution, p(w\r, 0°) is significantly skewed, in particular at large 
separations (EFWD; ZQSW; F94). The physical origin of the skewness and exponential 
shape of p(w) at large separations has until now remained unexplained. We provide the 
explanation below. 



2. The origin of the negative skewness 

Let v and 5 be the peculiar velocity of a galaxy and the mass density 
contrast at comoving position r 1; while v' and 5' the velocity of another galaxy at 
position r 2 at a certain fixed separation r = r 2 — r x . In our coordinate system 
r = {x,y,z} = r {sm8 cos ip, sin 9 sin ip, cos#}; the unit vector z points along the observer's 
line of sight, while v z and v' z stand for the line of sight velocity components. The 
probability that the four considered random fields reach values 3? = {8, 8',v, v'} in the 
range ddt = d5 d8' dw dv' is g(dt) ddt, and we will use brackets to denote ensemble averaging, 
(...)= J ... g d$l. As usual, expectation values (...) are assumed to be equal to spatial 
averages. The latter, however, should not be confused with number-weighted averages 
carried over galaxy positions (cf. LSS and F95). The iV-th moment of the relative velocity, 
weighted by the density of pairs of galaxies, is given by 

({v' z ~v z ) N {l + 8 g ){l + 5' g )) 

mN ~ <(i + W + ^)> ' () 

where 5 g = (5p/p) g is the contrast in the number density of galaxies, which may differ from 
the spatial fluctuations in the mass distribution. Here, we ignore this potential difficulty 
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and implicitly assume S g (r) = S(r). A central moment of order N is given by 

A* = ((v' z -v z - mi ) N (l + 5)(l + 5') ) [1 + e(r)]- 1 (2) 

To calculate the first few moments, we shall expand the random fields in perturbative series, 
v = v (1) +v (2) + . . ., 5 = (5 (1) +<5 (2) + . . ., where each superscript describes the perturbative order 
(5 {2) = 0[5 {1) ] 2 , etc.). We assume that all linear order terms are described by a joint-normal 
probability distribution. To lowest non-vanishing order, rrii(r,9) = v 2 i{r) cos# + 0(£ 2 ); 
fi a (r,d) = cr 21 2 (r,9) + 0(£ 2 ); /i 3 (r,9) = 6 ( (v™ - 2v^ ') ) + 0(£ 3 ), where 
£(r) = (5 W 5 {1) ') is the linear correlation function, while the mean streaming velocity, v 21 , 
and the pairwise velocity dispersion, <r 21 2 , are given by 

v 21 (r) = -2F/(0) f i{s){s/r) 2 d Sl (3) 

n(0) - n(r) cos 2 - E(r) sin 2 6] , (4) 

where f(£i) ~ fi ' 6 , while II and S are the radial and transverse components of the velocity 
correlation tensor, related to £(r) by equations (21.72)-(21.75) in PPC (see also |G6rski] 
1988| and |Groth et al. 1989|) . Note that the leading terms in the expansions for the first 



two moments come from linear perturbation theory. The third moment is different. In 
the early universe, linear perturbation theory is sufficient, and to first order, /i 3 = 0, 
in agreement with the assumed Gaussian initial conditions, symmetric about 5 = 0. 
However, gravitational instability breaks the initial symmetry (LSS; |Juszkiewicz et al. 1993 



Bcrnardeau et al. 1995| , hereafter BJDB). To calculate /i 3 , we need the second-order term 



for the velocity field. It can be obtained by inverting the expression for V • v (2) , derived by 
BJDB. For a curl-free flow, the resulting skewness is 

fi 3 (r, 9) = o 21 2 (S v + S A ) cos 9 , (5) 

where the first term, 

S v = 3v 21 {r) (cos 2 9 + C sin 2 9); C(Q) w f Q^ 21 , (6) 
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is induced by the mean streaming, while the second, 

S A = 6 [£(r) - n(r)] (1 + cos 2 6) (1 - C)/(Hfr), (7) 

comes from the anisotropy of the relative velocity dispersion tensor. When r _L z (both 
galaxies are in the plane of the sky), as well as for zero separation, the skewness vanishes: 
/i 3 (r, 90°) = fi 3 (0,6) = 0, in agreement with symmetry requirements. Such requirements do 
not apply to the radial component of the distribution (r || z) when the mean infall velocity 
is different from zero. We will now show that if £ remains positive for separations < r, 
then fi 3 (r, 0°) < 0. First, note that a positive £ implies v 2 i{r) < and S(r) > U(r) (see 
eq. [§ above and eqs. [21.72]-[21.75] in PPC). Hence, Sy(r, 0°) < while S A (r,0°) > 0: the 
effect of the negative infall velocity is counterbalanced by the anisotropy of the velocity 
correlation tensor (in agreement with N-body results; see p. 938 in F94). However, the two 
effects do not cancel out and for all 'reasonable' models the term, induced by the mean 
streaming dominates. For a power-law correlation function, £(r) oc r -7 , the ratio of the two 
terms is S A (r, 0°)/S v (r, 0°) = -8/7(5 - 7) « -0.35 for 7 = 1.7. 

In order to test our perturbative predictions against fully nonlinear N-body experiments, 
we used data, generated from the simulations described in Frenk et al. (1990; appropriate 
codes were kindly provided to us by Marc Davis). The simulations are of a standard CDM 
model with Q — 1, h — 0.5, and cr 8 = 0.6 (here cig is the rms density contrast in a 8 
h^ 1 Mpc sphere). These simulations contain N= 64 3 particles in a L = 180 h~ l Mpc box. 
The first three moments at separation r = 10.4 /i _1 Mpc, and 9 = 0°, determined directly 
from the simulations, are (t> 21 , cr 21 , ij. 3 1 ^) = (—190, 440, —400) km s -1 , while equations (|3|), 
(|j) and (|Sp give (—200,430,-430) km s _1 , respectively. Comparisons with numerical 
experiments of higher resolution (N = 256 3 , [Springcl et al. 1998 ) suggest that the accuracy 



of the N = 64 3 results, quoted above, is ~ 20%. We conclude that the perturbative 
predictions are in excellent agreement with the N-body experiments. 
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3. The origin of the exponential tails 



Eulerian perturbation theory, truncated at second order, can be used to write p(w\r, 9) 
as a marginal probability, obtained by integrating a 14- dimensional Gaussian distribution. 
Here we will not do that, however. Instead, we will trade accuracy for simplicity of 
calculations and replace rigorous perturbation theory with a following Ansatz. Suppose 
that the relation between the pairwise velocity and the random vector 3? = {5, 5', v, v'} is 
given by the mapping 



where u = v^' — and A = <5 (1)/ + 5 (1) . The first three moments of this new variable 
can be readily expressed in terms of v 2l and a 21 . To lowest non- vanishing order, we obtain 
(w) = v 21 (r) cos6*; ((w — (w)) 2 ) = cr 21 2 (r, and 



Clearly, the transverse (6 = 90°) components of the above moments agree with those 
obtained from rigorous second-order Eulerian theory in § |2] above. Hence, our Ansatz 
provides an acceptable approximation of the second-order prediction for p(w\r, 90°). For the 
radial part, we recover the true values of the first two moments only. The third moment, 
obtained from the approximate distribution is an overestimate of the true |/x 3 (r, 0°)|. 
However, the approximate expression @ correctly reproduces the negative sign of /z 3 , the 
scaling with the cosmological density parameter, /i 3 oc f2 L8 , as well as the scaling with the 
two lower moments, introduced by the dominant, infall term oc v 21 a 21 2 . 

According to equation (|]), w is the sum of two variables, u and w = uA. The 
probability distribution for u is 



w = u(l + 5 <1) ')(l + 5 <1) ) 



u + uA + 0(u 3 ) , 




((w-(w)f) = 6a 21 2 (r,9)v 21 cosfl . 



(9) 




(10) 
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The probability distribution for w = uA is readily obtained by integrating the expression 

+00 

p„(m\r,8)= / — tP u a(m, — J , (11) 
J \u\ \ uj 

—00 

where p u/A (u,zu/u) is a standard, joint-normal distribution for u and A (eg. F95), with 
w/u substituted for A. This integral gives 

p w {w\r, 0) = — K {(3\w\) , (12) 
na 21 

where K is the usual modified Bessel function, a = l/a A y/l — k 2 , (3 = l/a 21 a A (l — k 2 ), 
k = v 21 cos6/a A a 21 , and a\ = (A 2 ) = 2 [£(0) +£(r)]. There are several remarks worth 
noting about the properties of the random variable w. Near the origin, its distribution 
has a cusp; for small values of the argument, K (|a;|) = — [ln(|x|/2) + 0.577] + 0(x 2 ). 



For large values of the argument, Ko(|sc|) = yrc/2\x\ e - ^' [1 + C*(l/|a;|)]: this is the 
exponential behavior in the wings, typical for products of Gaussian fields (eg. the x 2 
distribution; see also |Scherrer 1992| or [Holzer and Siggia 1993| ). Finally, note that the 
exponential in equation QT2]) is not symmetric about w — 0, and its skewness is introduced 
by cross-correlation between the velocity and density (i.e., by the infall). The distribution 
of w = w + u is qualitatively similar to p ro ; it can be obtained from an expression, similar 
to eq. (|TTD with the original integrand, replaced by |m| _1 PmA [u, (w/u) — 1]. The resulting 
integral can be rewritten as 



p(w\r,6) = — exp I -— 
J a 2cr 

n v 



W(w,v 21 ,a 21 ) , (13) 



with 



na 21 



i w a 
a [ npa 

a 



(14) 



One can approximate the above integral via the method of steepest descent. This yields a 
convenient analytic formula for the probability distribution, 



cosh ^\U\/Ka A j \u\ a 2 \ 

p{w\r,9)fa ; — exp<^ — - — > (15) 

27rcr 21 a A \w\ { n z ) 
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where U = w/cr 21 a A , K = 1 + sgn(w)n, sgn(u>) = +1 for w > 0, and sgn(w) = — 1 for w < 0. 
The characteristic function is given by the Fourier transform of equation (pf), 

(e«*) = [0(t)]- 1 /2 eX p(- ( r 21 ¥/20(t)) , (16) 

where <f>(t) = 1 - 2ina 21 a A t + er*(l - K 2 )a 21 H 2 . In the limit ct a -> 0, (e iwt ) -> exp(-<x 21 2 t 2 / 2 ), 
and we simply recover the Gaussian distribution p u (w) (eq. [TO]) , in agreement with F95. This 
is obvious when we recall that non-Gaussian behavior of the distribution arises from the 
quadratic nonlinearities introduced by the number weighting; these quadratic terms become 
small compared to the volume weighted (and Gaussian distributed) velocity difference as 
the amplitude of the fluctuations decreases. In the opposite limit, when a A is increased, 
the distribution rapidly develops a central cusp and exponential wings. It is interesting to 
compare our results, valid in the transition zone between the linear and nonlinear regime, on 
separations r £ 5 h Mpc, with the analysis of Diaferio & Geller (1996) and Sheth (1996), 
restricted to much smaller scales, where a A ^> 1 and the perturbative approach breaks 
down. Note that the integral in equation ( fL3f) is a weighted sum over Gaussians having 
a range of dispersions. The weighting factor W is determined by the velocity correlation 
tensor and the velocity-density cross correlation, v 21 . This expression is similar to equation 
(5) of Sheth (1996), valid in the strongly nonlinear regime, where W is related to the 
Press-Schechter multiplicity function. The outcome of summing Gaussian distributions in 
both cases is an exponential distribution. A significant difference is that in the strongly 
nonlinear regime all sources of the velocity skewness vanish: the limit r — > 0, a A 3> 1 
corresponds to virialized cluster centers; there is no infall (v 21 = 0); the velocity dispersion 
is isotropic (II — S = 0), and /x 3 = by symmetry. 

We will now compare the predicted velocity distribution with direct measurements from 
simulations. Since the N-body results we have at our disposal assume a CDM spectrum, we 
need to introduce a shortwave cutoff in the initial power spectrum; otherwise £(0) becomes 
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infinite. We use a Gaussian filter of width R s and multiply the linear power spectrum, 
P(k) = J ^(r)exp(ik- r)d 3 r by exp(— k 2 R 2 s ). The resulting £(r) is finite at r = and 
remains flat for r £j R s - Existence of such a 'core radius' in any realistic £(r) is necessary 
anyway since galaxies are not point-like objects and have finite sizes. R s can also be related 
to the effective dynamical resolution of the simulations; we postpone the discussion of this 
problem to a later paper (Springel et al. 1998). Finally, the small-scale cutoff can be useful 
as a makeshift solution for reducing the \/i 3 \, overestimated by our Ansatz, and bringing this 
moment closer to the value predicted by the rigorous perturbative calculation in §2 above. 
We tried several filtering widths and found that R s = 3 /i _1 Mpc provides a reasonable fit 
to N-body simulations. In Figure 1 we compare the probability distribution, p(w\r,0°)dw, 
calculated from equation (13|), with direct measurements from N-body simulations. The 
upper panel shows the results of measurements from the simulations of Frenk et al. (1990; N 
Ri 2.6 x 10 5 ), while the lower panel - those from Zurek et al. (1994; N m 1.7 x 10 7 ), obtained 
by fitting their Fig. (7d) by a double exponential. The separation is respectively r = 10.4 
/i _1 Mpc, and r = 10.5 /i _1 Mpc. The sizes of the velocity bins are dw = 72 and 100 km s _1 
for the upper and lower panel, respectively. The assumed P(k) in both cases is the standard 
CDM spectrum, described in §2 (the only difference is that here we use P{k) exp(—k 2 R 2 ) 
with R s = 3 /i _1 Mpc, while in §2 R s = 0). Clearly, the agreement between the perturbative 
predictions and N-body measurements improves when the resolution of the simulations is 
improved. 

A possible alternative to the Ansatz, adopted in this section, is to derive p(w) from 
the |Zerdovich ( 1970 )| approximation (hereafter ZA). Like our Ansatz, this approach makes 



calculations simpler than the rigorous treatment in §2 above, and the simplicity here, too, 
is bought at the expense of accuracy in estimating /i 3 : at second order, the ZA breaks the 
momentum conservation flJuszkiewicz et al. 1993| ; BJDB) by implying C = in eqs. @ 



and (|7|). As a result, the ZA underestimates |/i 3 | by ~ 50%. At first order, however, 
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the ZA agrees with the rigorous Eulerian perturbation theory, so its predictions for v 31 
and <7 21 must be identical with ours. The ZA was recently used by |Seto and Yokoyama 



1998)| . Qualitatively, their p(w) is similar to ours. However, at the quantitative level we 



disagree because their method underestimates v 21 by at least an order of magnitude. As a 
consequence, Seto and Yokoyama had to readjust their predicted p(w) 'by hand' to achieve 
agreement with simulations. Their results seem puzzling given the properties of the ZA, 
discussed above. 

Another alternative approach one might consider, is to expand p(w) in orthogonal 
polynomials (eg. |Juszkiewicz et al. 1995| ; |hifshitz and Pitaevskii 198 1| , p. 31). At the end, 



only direct applications to future redshift surveys, like the SDSS or 2dF, will decide which 
of the discussed physical models of p{w) provides the optimal combination of simplicity and 
accuracy. 
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FIGURE CAPTION 
Figure 1 

Comparison of the probability distribution, dP = p(w\r,0°)dw, derived from 
equation ( |TB| ) (shown as curves) with direct determinations from N-body experiments 
(shown as histograms). The number of particles in the two sets of simulations was N = 
2.6 x 10 5 (upper panel) and N = 1.7 x 10 7 (lower panel) . Note that an increase in resolution 
brings the N-body results closer to the perturbative predictions. 
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